#include "bc_freeOutflow_velocity.hpp"
#include "meshTypes.hpp"
#include "field.hpp"
#include "boundaryPatch.hpp"
#include "simulationTime.hpp"
#include "exceptions.hpp"

namespace hms {

template<typename MeshType>
FreeOutflow_velocity<MeshType>::FreeOutflow_velocity(
	 Field<MeshType>& v, const BoundaryPatch& patch
) : BaseType(Type::freeOutflow_velocity, v, patch)
{}

template<typename MeshType>
void FreeOutflow_velocity<MeshType>::update() const {
	Field<MeshType>& v { field() };
	const BoundaryPatch& p { patch() };
	/* for ingoing flux, this BC behaves like a solid wall: both
	 * velocity components are inverted.
	 * Otherwise, both components are copied from the domain values
	 * So we form the projection of vDomain onto the outward normal.
	 * Its sign is positive for outflow and negative for inflow.
	 * Multiplying this sign with vDomain gives vGhost.
	 * */
	/* rectangular meshes */
	if constexpr (
		std::is_same_v<MeshType, UniMesh> ||
		std::is_same_v<MeshType, RectMesh>
	){
		using BL = BoundaryLocation;
		BL location { whereIs( p, v.mesh().size() ) };

		/* for the top and right boundary, a positive value of the velocity
		 * component orthogonal to the boundary edges signifies outflow,
		 * whereas the same value signifies inflow at the bottom and left. */
		int sign     {location == BL::bottom || location == BL::left ? -1 : 1};
		int orthoRow {location == BL::bottom || location == BL::top  ?  1 : 0};

		v.ghostValues(p) = sign *
			v.domainValues(p).row(orthoRow).sign().template replicate<2,1>() *
			v.domainValues(p);
	} else {
		v.ghostValues(p) = 
			( // projection is colwise sum of coeff-wise product
				v.domainValues(p) *
				v.mesh().edges().outwardNormals(p).array()
			).colwise().sum()
				.sign().template replicate<2,1>() *
			v.domainValues(p);
	}
}

#undef CHEAD

#define INSTANTIATE_FOR_ALL_MESHTYPES(MeshType) \
template class FreeOutflow_velocity<MeshType>;
#include "meshInstantiator.hpp"


} //end of namespace
